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In high energy heavy ion collisions of RHIC and LHC, a strongly interacting 
quark gluon plasma (sQGP) is created. This medium undergoes a hydrody¬ 
namic evolution, before it freezes out to form a hadronic matter. The initial 
state of the sQGP is determined by the initial distribution of the participat¬ 
ing nucleons and their interactions. Due to the finite number of nucleons, the 
initial distribution fluctuates on an event-by-event basis. The transverse plane 
anisotropy of the initial state can be translated into a series of anisotropy 
coefficients or eccentricities: second, third, fourth-order anisotropy etc. These 
anisotropies then evolve in time, and result in measurable momentum-space 
anisotropies, to be measured with respect to their respective symmetry planes. 
In this paper we investigate the time evolution of the anisotropies. With a nu¬ 
merical hydrodynamic code, we analyze how the speed of sound and viscosity 
influence this evolution. 


1. Hydrodynamics in ultra-relativistic heavy ion collisions 

The medium created in ultra-relativistic heavy ion collisions is success¬ 
fully modeled by perfect fluid hydrodynamics, i.e. the observables of soft 
hadrons, photons and leptons can described by hydro models. Exact solu¬ 
tions are particularly interesting, since one gains an analytic insight on the 
connection between the initial and final state of the hydro evolution. 

These solutions fulfill a simple set of laws, governing many systems on 
many scales. Hydrodynamics only assumes the local conservation of energy 
and momentum, and, depending on the circumstances, of entropy or some 
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conserved charges. The equations of non-relativistic hydrodynamics are: 

+ Vpt; = 0 (1) 

p{^ + (^’^)^) = -^P + + (^ + f) (2) 

dc 

^ + Vev = —pVv + V(av) (3) 

at 

where p is matter density, v is the velocity field, s is the (internal) energy 
density, p is pressure, a is the viscous stress tensor of the homogeneous 
and isotropic medium, p is the dynamic viscosity coefficient, A the bulk 
viscosity coefficient, while t represents time and V = (dx^dy^dz) the spa¬ 
tial derivative vector. If we however would like to utilize hydrodynamics 
in a relativistically expanding medium, we have to use the equations of 
relativistic hydrodynamics. In case of perfect hydro, these take the form of 

(e + p)—(4) 

To close this set of equations, we need to specify a relation between energy 
density e and pressure p. This we obtain using the below Equation of State: 

5 = k{T)p, (5) 

where k = c~‘^ is the inverse square of speed of sound. It’s temperature 
dependence can be taken from lattice QCD calculations, but often n is 
simply kept constant. 


2. Initial conditions and anisotropies 


There are many solutions that fulfill the above equations, we however need 
to start from an initial state that corresponds to the energy distribution 
deposited by the nucleons participating in the collision. In nucleus-nucleus 
collisions the initially formed medium can be spatially approximated by an 
ellipsoid, which can be described by the contour levels of a scale variable s 


c =_L ^_^_ 

Z2 


(6) 


where X, T, Z are the axes of the ellipsoid described by the equation s = 
1. Multiple exact solutions of hydrodynamics were foundllll that resemble 
this symmetry, i.e. where the thermodynamical quantities are constant on 
s =const. curves. It is important to see however, that nuclei contain a finite 
number of nucleons, and thus the participants of a nucleus-nucleus collision 
form an event-by-event fluctuating shape (see Fig. [^. 











February 16, 2016 2:23 


WSPC - Proceedings Trim Size: 9in x 6in 


bagolycsanad_numhidro 


3 



Figure 1. 


Simulated participant nucleon distribution (figure from Ref0) 


Thus a simple ellipsoidal symmetry is not sufficient to describe the hy- 
drodynamically expanding medium, not to mention that it does not de¬ 
scribe higher order flow coefficients either. To overcome this difficulty, one 
can assume higher order anisotropies in the scaling variable, by redefining 
it as 

2 2 

s = ^ (1 + S 2 cos(2(/))) + ^ (7) 

with r and (p being the polar coordinates in the transverse plane, while R 
and 82 are in a direct relation to X and Y. This can then be generalized to 




n=2,3,. 


e„ cos mp\+^ 


(8) 


in the cos rup expression a \l/„ n-th order event plane can also be introduced. 
Example distributions generated with this scaling variable are shown in 
Fig.H Solutions utilizing scaling variables like the above s were presented 
in Ref. SI However, these solutions, as well as the ones in Ref. SI are valid 
only in case of a 3D Hubble-flow, which corresponds to the lack of pressure 
gradients. This also means that the anisotropies will remain constant in 
time, which may not be the case throughout the whole evolution. In our 
analysis, we thus apply numerical calculations to include pressure gradients 
in the initial conditions, and analyze the time evolution of the anisotropies. 

Thus we will start the hydrodynamic evolution from the above described 
set of initial conditions. We would like to extract the anisotropy coefficients 
at later stages of the evolution, for which we need a definition. We in¬ 
troduced the anisotropies of density, energy density and velocity fields as 
follows: 


— 


(cos(n<^))p/£/w, 


(9) 
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Figure 2. Initial distribution with the outlined s scaling variable, for different Sn values. 


where w = exp {—v‘^ ~ defined to calculate the asymmetry of velocity 

field. To understand the relation of the asymmetry parameters introduced 
byi to the scale variable anisotripy coefficients, we checked their 

relation. As a first order approximation (in Sn), we obtained 

e, = -Sn/2 (10) 


however, in a more accurate approximation the result is: 


ei p 

(^2 + £4)^3 _ 

2 + Hn 

+ + 0 ( 4 ) 

(11) 

62 - 

—62 + ^ 2^4 

n 

(12) 

63 - 

“^3 

-y + ^£3 d + 

n 

(13) 

2 + Fn 

64 P 

c- _L lc -2 
c 4 “T 2^2 

-y + ^ - 4^4 XI d + C>(e^) 

n 

(14) 

2 + Fn 


this means that even if there is no Si coefficient, we will still get an ei. 
It is also important to note that ^2 7 ^ 0 and = 0 will result in a 

nonzero 64 = 5 ^/ 4 . See an example with multiple nonzero Sn coefficients in 
table where we give an example how a set of values result in a different 
set of Cn values. The actual results can be accurately approximated with 
eqs. @-([T4]). 


3. Numerical method 

Let us now discuss how we evolve the previously defined initial conditions 
in time. In order to do that, we first transform the equations of (1+2 di¬ 
mensional) hydrodynamics to a convective form as usual: 


StQ + a,F(Q)+a^G(Q) = 0 


( 15 ) 
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Table 1. Example for the — Sn relation. 


n 

£n 

direct calculation 

en 

approximation 

first order estimate 

1 

0 

0.00278 

0.00298 

0 

2 

0.1 

-0.04927 

-0.04869 

-0.05 

3 

0.05 

-0.02533 

-0.02484 

-0.025 

4 

0.02 

-0.00769 

-0.00745 

-0.01 


where Q is calculated from the hydrodynamical fields, while F and G are 
functions of these, transformed from the original hydro equations. The finite 
volume method discussed below and introduced in Ref.l^ works for many 
different systems of partial differential equations (PDE’s), one just needs 
to determine the given functions to transform the PDE’s to the above form. 
This can be done in a straightforward manner, see for example RefsP^ 

The basis of finite volume methods is to discretize the Q vector, by 
defining it on a space-time grid, and obtaining the discrete values by av¬ 
eraging in each cell. Then, we have to evaluate fluxes F and G between 
grid points (these are then called intercell fluxes). This is difficult, since the 
fluxes depend on the values of Q, but these are defined on the grid points. 
Our method, outlined in Ref.l^ and used for exampe in Ref,l^ is based on an 
accurate multi-stage (MUSTA) predictor-corrector estimation of the inter¬ 
cell fluxes. We discuss this method briefly in 1+1 dimensions (and with a 
single-component Q only) - it is straightforward to generalize it to multiple 
spatial dimensions via operator splitting (we used the Lie type of splitting 
method here). In case of viscous hydrodynamics, we also used operator 
splitting to separate ideal and viscous fluxes. 

To estimate the intercell fluxes between cells i and i + 1, we start from 
the Q2 and values, where n represents the current time-step. Let us 
call these values indicating that these represent a “zeroth order” 

estimate of the intercell values. We will then make a series of estimates in 
the given cell, the Ith estimate being indicated with The flux in the 

cell centers is then = F We define an intermediate Q value 

and flux F as follows: 




+ Qfli 


2 Aa; L 



F«=f(q«J (16) 


Then our intercell flux estimate is 


i+ o 


p(l) I or(0 




Ax 

At 


(e+-+) 


(17) 
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This is the so called FORCE flux approximation. The essence of our the 
applied MUSTA algorithm is that now we make a new prediction, 
and with that we repeat the steps and we get a better flux approximation, 
and with that we can make better prediction, and again better flux approx¬ 
imation. To make a new prediction simply we use the equations we want 
to solve: 


q(*+1) ^ qW . 

= Qfh 


At 
Ax . 
At r 
Ax 


i{i) 

^i+1 




771(0 


(18) 

(19) 


We may then stop at a given I value, by setting i With this, 

we simply calculate the next time-step as 


<2o> 

where F^_i can also be calculated by the procedure outlined above, if we 

* 2 

start from cells i — 1 and i. 

In order to complete the method, let us mention that we used the initial 
conditions as described above, we used non-reflective boundary conditions, 
and the timestep was determined adaptively via the Courant-Friedrichs- 
Lewy condition. We tested our method with the usual Sod shock tube and 
known hydrodynamic solutions as well. 


4. Results 

With the above outlined method, we solved the equations of nonrelativis- 
tic hydrodynamics with viscosity, as well as the equations of relativistic 
perfect hydrodynamics. Our goal was to investigate how various speed of 
sound values and various viscosity values affect the time evolution of the 
asymmetries. 

Let us first discuss the effects of viscosity (in nonrelativistic hydro). As 
shown in Fig. [^viscosity makes the disappearance of pressure anisotropies 
slower, and the same is in the number density. In other words, in case 
of nonzero viscosity the pressure asymmetries remain in the system for a 
longer time. This may be due to the fact that viscosity makes the flow 
slower, as illustrated in the first two plots of Fig. As shown in the bot¬ 
tom plot of Fig. and the third and fourth plots in Fig. the effect of 
viscosity on the flow field is however the opposite: viscosity makes the flow 
anisotropies disappear faster. The reason of this may be that the viscous 
force contains the second spatial derivatives, thus the fluid elements with 
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larger asymmetries feel a stronger viscous force, which means that viscosity 
“washes out” the anisotropies. 

We also investigated the effects of a change in the speed of sound. As 
shown in Fig. the reduction of speed of sound makes the anisotropies 
disappear slower. The reason for this is that pressure waves travel with the 
speed of sound, so the equalization of pressure is slower in case of smaller 
speed of sound values. This is confirmed by relativistic calculations as well. 
As shown in the top panels of Fig. in case of a larger k (i.e. a smaller 
value) the anisotropies disappear slower, and also the cooling is slower. 
This also means that of course the system will freeze out later. The time 
evolution of the flow field is shown in the bottom panels of Fig. 



t [fm/c] 



Figure 3. Time evolution of pressure anisotropies is shown in the left panel, while 
anisotropies in the flow field are shown in the right panel. 


5. Summary 

In this paper we analyzed the time evolution of the asymmetries of the 
QCD matter created in high energy heavy ion collisions. We utilized nu¬ 
merical hydrodynamics instead of exact solutions to investigate situations 
not describable by the presently known hydrodynamic solutions. We how¬ 
ever did not start from the most realistic initial conditions, but started 
instead from an initial condition that is very similar to one described by 
known analytic solutions - except in pressure, where we used a pressure 
profile similar to the density profile given in these models. We analyzed 
both non-relativistic and relativistic hydrodynamics, and arrived at similar 
conclusions. It turns out that a stiffer equation of state, i.e. a larger speed of 
sound makes the pressure anisotropies disappear faster. The appearance of 
viscosity also makes the flow anisotropies disappear faster, however, pres- 
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Figure 4. Time evolution of the energy density (first two rows) and of the flow fleld 
(last two rows). The first and third rows show the viscosity free case, the second and 
fourth rows show the viscous case with fi = lOMeV • fm. 


sure anisotropies remain in the system longer, due to the slower flow of the 
viscous system. 
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Figure 5. Time evolution of pressure anisotropies is shown in the left panel, while 
anisotropies in the flow field are shown in the right panel. 


ful discussions. This manuscript was supported by the OTKA grant NK 
101438. 


Bibliography 

1. T. Csorgo, L. P. Csernai, Y. Hama and T. Kodama, Heavy Ion Phys. A21, 73 
(2004). 

2. T. Csorgo, S. V. Akkelin, Y. Hama, B. Lukacs and Y. M. Sinyukov, Phys. Rev. 
C67, p. 034904 (2003). 

3. C. Loizides, J. Nagle and P. Steinberg (2014). 

4. M. Csanad and A. Szabo, Phys.Rev. C90, p. 054911 (2014). 

5. E. F. Toro and V. A. Titarev, J. Comp. Phys. 216, 403 (2006). 

6. H. Nishikawa and K. Kitamura, J. Comp. Phys. 227, 2560 (2008). 

7. M. Takamoto and S.-i. Inutsuka, J.Comput.Phys. 230, 7002 (2011). 

8. I. Karpenko, P. Huovinen and M. Bleicher, Comput.Phys.Commun. 185, 3016 
(2014). 


























February 16, 2016 2:23 


WSPC - Proceedings Trim Size: 9in x 6in bagolycsanad_numhidro 


10 




t [fm/c] 


t [fm/c] 



Figure 6. The time evolution of velocity field. The top row shows k = 2 case, the 
bottom row the n = A case. Time evolution of pressure anisotropies is shown in the left 
panel, while anisotropies in the fiow field are shown in the right panel. 












































































































































